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^ (54) Title: DETECTION OF FEATURES IN IMAGES 

(57) Abstract: An image processing technique which identifies pixels in images which are associated with features having a selected 
shape, such as but not exclusively step edge, roof, ridge or valley. The shape of the intensity profile in the image is compared in 
an intensity independent way with a shape model to select those pixels which satisfy the shape model and are thus associated with 
the feamre of interest. This comparison is achieved by examining the phase and amplitude of a spectral decomposinon of parts of 
the image profile in the spatial or spatio temporal frequency domain. This decomposition can be achieved using quadrature wavelet 
pairs such as log-Gabor wavelets. The difference between the odd and even components, known as the feature asymmetry, gives an 

indication of the shape of the feature. The analysis may be extended to the time domain by looking at the shape of the image profile 

across a time sequence of images, which gives an indication of the velocity of a moving feature. Pixels identified as belonging to 
a feature of the right shape are labelled with the value of feature asynunetiy, the local amplitude, feature orientation and feature 
velocity, and this information can be used to improve the tracking of detected features through a sequence of images. 
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DETECTION OF FEATURES TN IMAGES 

The present invention relates to a method and apparatus for processing images to 
detect particular features in the image. It is particularly concerned with the detection 
of features in noisy 2D, 2D + T, 3D and 3D + T images, and is thus particularly useful 
in medical imaging, for instance using ultrasound. 

There are a variety of techniques for processing images to .enhance them, particularly 
to reduce the amount of noise and/or to detect features of interest in the image. Such 
techniques are particularly useful in areas where images are noisy, for inistance in 
medical imaging. For example, ultrasound images include a considerable amount of 
speckle which makes images difficxilt to inteq^ret Skilled operators are able to 
int^ret such images based on their experience and knowledge of the features of 
interest likely to be present. However, automating the process of detecting features of 
interest is made difficult by the presence of noise and imaging artefacts. Similar 
problems arise in other imaging techniques, such as magnetic resonance imaging. 

In general up to now techniques proposed for enhancing images have been based on 
examining the intensity of the image, for instance by smoothing the intensity to 
reduce noise or by setting a detection threshold based on amplitude. However, 
techniques based on thresholding are not entirely satisfactory because of the difficulty 
of setting a suitable global threshold. For instance, ultrasound images include 
attenuation artifacts caused by the changix^ orientation and reflectivity, of the tissue 
being imaged. Choosing one threshold v^ch is suitable for excluding noise, but 
including all features of interest, is difficult if not impossible. 

Other techniques have been proposed for analysing noisy images for detecting 
features which are based on correlating images over time or with movement of the 
imaging probe. For instance, ultrasound speckle decorrelates with movement of the 
probe and over time. A problem with some of these techniques, however, is that they 
tend to blur the features of interest in the image. Thus the user has to accept to trade- 
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off between removing noise and losing image quality. 

It is therefore an object of the invention to provide a technique for automatic detection 
of features of interest in an image which improves on the prior art techniques. 

5 

The present invention provides a method for detecting features of interest in an image 
based on the shape of the intensity profile of the feature, rather than its intensity. This 
is achieved by an intensity independent comparison of the intensity profile across the 
image with a shape model of features of interest To improve subsequent processing 
10 of the image, areas of the image which are detected as corresponding to the shape 
model are labelled with a description of their properties. 

In more detail, therefore, the present invention provides a method of detecting 
features-of-interest in an image, comprising the steps of: 
15 making an intensity independent comparison of the shape of regions of the 

intensity profile across tiie image with a shape model for predefined features of 
interest to define as features of interest areas of the image satisfying the shape model; 
and 

storing for each of said defined areas a label comprising the level of similarity 
20 to the shape model. 

The predefined features of interest may be positive and negative step edges, or roof or 
ridge shapes, or valleys or other features as desired depending on the application. 

25 The label may include a measure of the orientation of the feature in the image, and the 
comparison may also be performed on an intensity profile taken across successive 
firames of an image sequence to derive a measure of the velocity of image features. 
The label may then comprise a measure of the velocity. 

30 The comparison with the shape model is advantageously performed in the spatial or 
spatio temporal firequency domains, thus by decomposing the intensity profile into 
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spatial frequency components and examining the phase and amplitude of those 
components. In one embodiment this is achieved by convolving the intensity profile 
with a quadrature filter to derive the phase and amplitude of a pair of odd and even 
* components of the profile in the spatial frequency domain. The phase and amplitude 

5 of these components is characteristic of different shapes in the intensity profile. The 
difference between the odd and even components, which is a measure of the ^^eature 
asymmetry" is, for example, a maximum for a step edge. Thus by setting constraints 
on the value of the feature asymmetry it is possible to determine that the area of the 
image under consideration corresponds to the shape being sought. The values of 

10 feature asymmetry and local amplitude (which is based on the amplitude of the two 
components), are included in the label. 

Preferably the filters are quadrature wavelet filters, e.g. log-Gabor filters, so that the 
intensity profile is convolved with odd and even wavelets. The scale of the wavelet is 
1 5 chosen in accordance with the scale of the features to be detected, and to be much 
larger than thcscale of noise in the image. This means that the technique selects 
features according to their shape and scale, but regardless of the value of intensity. 

The filters may be oriented in different directions across the image to detect image 
20 features having different orientations. The label may then include a measure of the 

feature orientation, v^ch can be derived from the relative responses of the differentiy 
oriented filters. 

The image processing technique of the invention may be used as a precursor to 
25 tracking detected image features through a sequence of image frames. The provision 
. of the labels carrying information about the detected features is particularly usefiil in 
such a tracking process. For instance, rather than simply searching for the closest 
* image feature in two successive frames, it is possible also to compare the labels of the 
detected image features in the two frames and to define them as relating to the same 
30 image feature if the labels satisfy predefined search criteria. 



wo 02/43004 PCT/GBOl/05094 



-4- 

Such search criteria may include a coaditioa on the value of the feature asymmetry of 
the two features being detected and a condition on the orientation of the detected 
features (for instance that they are oriented similarly). 

5 Where the label includes a measure of the velocity of the feature, the search area from 
frame-to-frame may be defined in accordance with the velocity, and the search criteria 
may include a condition on the velocity. 

The invention is particularly applicable to ultrasoimd or magnetic resonance images, 
10 but also to x-ray and other imaging modalities. The embodiment described below 
relates to echocardiography, but the invention is applicable to ultrasound images in 
general, such as of different organs and tissues e.g. coronary arteries, liver, foetus etc, 
and to different ultrasoimd modalities such as the use of contrast and different signal 
processing techniques such as Doppler imaging and harmonic imaging. For instance 
15 in echocardiography, by setting the shape model ^propriately the invention is 

adaptable to the detection of any features corresponduig to a step change in intensity, 
such as the ventricular walls. 

The invention may be embodied in a system adapted to perforai the image processing 
20 method, and it may be embodied by a computer program comprising program code 
means for executing the method. Thus the invention further provides for a computer 
program storage mediimi cairyiag such a computer program, and also a computer 
system programmed to carry out the method. 

25 The invention will be further described by way of example with reference to the 
accompanying drawings in which:- 

Figures 1(a) and (b) illustrate Fourier decompositions of a step edge and a 
triangular function; 

30 

Figure 2 illustrates the values of the local phase signatures of spatial 
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firequency compoiients for diSbrent functions; 

Figures 3 (a) and (b) illustrate respectively even and odd two-octave log- 
Gabor wavelet functions; 

Figure 4 illustrates the spectral properties of the log-Gabor wavelets; 

Figure 5 (a) illustrates a typical intensity profile across an ultrasound image; 

Figure 5 (b) illustrates the local amplitude result of a convolution of an 
intensity profile with log-Gabor filters at a succession of scales; 

Figure 5 (c) illustrates the local phase result of a convolution of an intensity 
profile with log-Gabor filters at a succession of scales; 

Figure 6 schematically illustrates an image processing system according to an 
embodiment of the invention; 

Figures 7 (A) and (B) illustrate examples of the feature orientations in two and 
three dimensions used in the embodiment of Figure 6; 

Figure 8 illustrates the effect on feature detection of analysing the intensity 
profile using^log-Gabor wavelets of different scale (wavelengths); 

Figure 9 illustrates a process for tracking a detected feature through a 
sequence of images; 

Figure 10 schematically illustrates part of the tracidng process; and 

Figure 1 1 illustrates the relationship between local phase and feature 
asymmetry. 
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While traditional methods for detection of a step edge in an intensity profile rely on 
examining the derivatives of the intensity profile (i.e. its slope), to detect the sudden 
changes in intensity which mark the visible edge of features in the image, this 
embodiment of the invention examines a decomposition of the intensity profile in the 
5 spatial or spatio temporal firequency domains. For instance. Figure 1 illustrates 

Fourier decompositions of a step edge (Figure 1(a)) and a triangular function (Figure 
1(b)). In Figure 1 (a) it can be seen that three sinusoidal components i, ii and iii sum 
to make an approximation to a step edge iv. Further it will be noted that the phase of 
all three components is the same, zero, at the positive-going step edge. Further the 
1 0 phase value of all the components at the negative step edge is 1 80 degrees, jReference 
to Figure 1 (b) illustrates that for a triangular peak the phiase value of all of the 
Fourier components is 90 degrees at the peak of the triangle. 

The "local phase" is defined as the phase ofTset of each Fourier component describing 
15 a signal, measured at the event of interest. For features such as a step, the peak of a 
roof or the foot of a ramp, the local phase of all of the Fourier components converges 
to the same value. Further, the value is characteristic of the type of feature. The 
phase values for different types of feature is illxistrated in the diagram of Figure 2. It 
can be seen, therefore, that the phase of the spectral components of a signal gives an 
20 indication of the shape of the signal. It should be noted that the particular value of 
phase is arbitrary (depending on the position of the origin), what is important is that 
the different features have different, known values. 

While Fourier decomposition might be suitable for perfectly shaped features, the 
25 intensity profile across typical images is much more complex and noisy and a Fourier 
decomposition of such a complex waveform would be extremely complex, lostead, 
this embodiment of the invention looks in turn at successive small sections of the 
intensity profile across the image by making the spectral analysis using wavelets. This 
embodiment uses log-Gabor wavelets as illustrated in Figure 3 (which are in fact 
30 Gaussian modulated sinusoids) which have a non-zero central section and are zero 
elsewhere. Thus convolving them with the intensity profile eflfectively windows a 
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small region of the profile. Figure 3 (a) shows the even wavelet based on a Gaussian 
modulated cosine, and Figure 3 (b) the odd wavelet based on a Gaussian modulated 
sine function. The spectral properties (plotted against firequency and log-&equency) 
of the wavelets are shown in Figure 4. 

5 

Thus just as a Fourier analysis of a profile gives the amoimt and phase of each of the 
Fourier components, convolving the intensity profile with the log-Gabor wavelet pair 
extracts the amplitude and phase of the odd and even components for the region of 
interest. 

10 

The scale of the wavelets used is important Figure 5 (a) illustrates a typical intensity 
profile across an echocardiographic image, bne can imagine analysing this profile 
using very small scale versions of the wavelets shown in Figures 3 (a) and 3 (b). This 
would result in an analysis of the very small scale variations in intensity in the profile. 
1 5 On the other hand, if the scale of the wavelets used was much larger, such that, for 
example, only one wavelet fitted across the entire length of the profile, the effect of 
the small scale variations in intensity would not affect the result, and the analysis 
would only actually pick up the broad overall shape of the profile. Thus the effect of 
varying the scale of the wavelet is to look for different scale features in the profile. 

20 

The effect of varying the wavelength, or scale, of the wavelet is illiistrated by the 
scalograms of Figures 5 (b) and (c). These scalograms are constructed by taking the 
convolutions of a sequence of log-Gabor filters of increasing wavelength, and 
stacking vertically the output local phase or local amplitude responses of each filter, 
25 so that each forms a new row of the image. The scalograms on the left (Fig. 5B) are 
computed firom the local amplitude responses given by 

Aix)^^o\x) + e\x) 

where e(x) is the ou^\zt of the even log-Gabor wavelet, and o(x) is its odd 
30 counterpart. Similarly, the phase scalogram of Fig. 5B is computed Scorn the absolute 
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local phase values given by: 



0(x)= arctan 




On the left scalogram, the intensity values correspond to the output amplitudes of the 
convolutions. The phase values on the right have been mapped to intensities so that 
black corresponds to zero phase (positive steps) and white coiresponds to it radians 
(negative steps). Peaks, having 7c/2 phase signatures, are grey. 

The phase scalogram shows that at small scales many features correspondiag to these 
extreme values are found. However, as the filter scale increases, the number of 
features reduces to a few straight bands. Among these it is possible to identify the 
step edges (local phase converges to zero, thus black on the scalogram) corresponding 
to the epicardial (1) and the endocardial (2,3) step edges, as well as a pericardial peak 
(4). From the phase scalogram such as Fig. 5c, two observations can be made: 1) 
features such as step edges and peaks have stable phase values and 2) phase 
signatures remain localised over scales. Stability means that there is always a 
particular phase signature associated with each feature shape over scales (in the 
illustration the zero degree and it radian values corresponding to positive and negative 
steps respectively are shown). Localisation means that the local phase signature 
associated with a feature remains in the proximity of the event (i.e. at the same place 
on the X-axis) over a range of scales. Thus because localisation and stability hold, it 
is possible to identLfy features of interest from the phase output of a large scale 
convolution according to pre-defined phase signatures, and trace these back along 
straight paths to their location. Thus the system can use a large scale filter to remove 
noise and still identify feature points accurately, and by examining a phase scalogram 
such as Figure 5 (c)it is possible to set a suitable scale at which tiae phase signatures 
of features of interest are clear. 



In the present embodiment of the invention a single scale is sufBcient to detect echo 
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boundaries of interest corresponding to the heart, hi different ^plications, such as 
viewing different organs or using different modahties it may be necessary to combine 
the outpvrt of the filters at different scales to enhance the detection. If necessary filters 
may be run at different scales and the results of different scales combined by, for 
5 example, averaging feature asymmetry resxxlts and recording selected features 
if their shape description is consistent over a number of scales. 

In this embodiment of the present invention the images being processed are 
echocardiograms of the left ventricle obtained using an HP-Sonos 5500 scanner 
10 producing images of a size 300 x 250 pixels at a firame rate of 50 to 60 Hz. For these 
images scales of between 32 to 56 pixel are used with log-Gabor wavelets of two 
octaves. 

So far, the issue of boundary detection has been centred on the ID problem and the 
1 5 analysis of a ID profile. However, the present invention is also designed to detect 
and find the orientation of features (steps, valleys, ridges) in 2D image and 3D 
volume data (including 2D+T i.e. sequence data). 

To do this a filter bank of oriented filters is designed to sample the space of the data: 
20 i.e. a plane for a 2D image or a volume for 3D or 2DH-T data. This embodiment of the 
invention has been applied to echocardiograptdc sequences, and in this case a 2D+T 
yolxmae is created by stacking the successive image fi-ames behind each other creating 
a 3D volume (such that the two spatial dimensions follow the x and y axes and the 
time dimension follows the z axis). A volume for three spatial dimensions (3D) is 
25 created by stacking images taken at different depths. 

Sampling of multidimensional data and detection of feature orientation require: 

1 . the extension of log-Gabor filters to higher dimensions (2D and 3D), 

2. the choice of optimal orientations for the filters in the 2D and 3D spaces, and 
30 3. an interpolation method that allows estimation of the orientation of the 

features detected firom a small* number of filter responses. 
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Extending a ID filter function to higher dimensions is called 'spreading'. In this 
embodiment the log-Gabor filters are spread to higher dimensions by multiplying 
with a Gaussian fiinction in the spectral domain. For example in 2D 



G(a),,^)= exp- 



■ + 



Here Of> defines the extent of the spreading function as a scaling s of the separation 
between filters cr^= s x A ^ . This spreading concept can, of course, be extended to 
3D. 



The number of filters and their particular orientations have to be defined such that all 
possible features of interest are detected and their direction estimated reliably. In the 
. absence of prior knowledge on the expected orientation of the features, this requires 
evenly sampling the data space with sufficient filters, so that an interpolation 
procedure can estimate the orientation of the feature accurately. In this embodiment 
three evenly spaced filters in 2D and six evenly spaced filters in 2D+T provide good 
accuracy while minimising the number of filtering operations required. These are 
illustrated in Figures 7(A) and 7(B) respectively. 

Various interpolation algorithms can be used; in this embodiment, the filter responses 
are interpolated by using an elliptic (2D) or ellipsoidal (3D) fit to the filter responses. 
This approxin:iation provides estimates of feature orientation within the accuracy 
required by the tracking algorithm. 

Estimation of feature orientation in a 2D+T volume is equivalent to finding the 
orientation and direction of movement of a feature. If one imagines a horizontal Une 
moving upwards over a sequence of images, then it will appear in the frame stack as a 
plane with a given inclination to the horizontal proportional to the speed of the line. 
A fast line will become a steep plane; a slower line would be less steep. Estimating 
this component of feature orientation is equivalent to estimating feature velocity. 
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This can be used to anticipate the motion of anatomical boundaries and is of great 
benefit for subsequent tracidng stages. 

This embodiment of the invention is particularly applicable to the detection of 
anatomical boundaries such as heart chamber boundaries and in particlar ventricular 
boundaries in the echocardiographic image of the left venticle. These are step edges 
in intensity. In this embodiment such step edges are detected by deriving the "feature 
asirmmetry"' which is a measure of the difference between the odd and even filter 
outputs at each pouit m the image and is a function of the phase signature. Thus 
rather than lookii:]^ at the value of the local phase as defined above, it is the feature 
asymmetry which is derived and this can be defined as:- 

where o(x) and e(x) are the odd and even filter outputs at point x in the 2D image or 
2D+T volume. 

The feature asymmetry is directly related to the local phase in the following way:- 

FA(x) = |siii(<E)(x))| - |cos(<D(x))| 
where O is the local phase. 

This fionction is plotted in Figure 11. It can be seen that the feature asymmetry varies 
almost linearly with local phase ranging fiom -1 for peaks and valleys to' +1 for 
positive and negative steps. (It should be noted that the phase values for Fig. 1 1 and 
the embodiment below are offset by 90** firom those shown in Figs. 1 and 2. If the 
zero phase point had been defined as in Figures 1 and 2, then Fig. 1 1 would simply be 
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inverted, with FA = -1 for positive and negative steps (at 0° and 1 80**) and +1 for 
peaks and valleys (at 90° and 270°).) 

With this embodiment of the present invention, therefore, it is the feature asymmetry 
which is calculated and thresholded to detect positive and negative steps. For 
instance, a threshold of 0,9 for the feature asymmetry corresponds to a local phase 
range of ± 5 degrees around ± 90 degrees and will thus detect positive and negative 
steps. Theexact value ofthe threshold affects the width of the detected features. The 
threshold is dimensionless and thus can be applied to any image to detect the edges 
that fall within the scale of the filters. Noise rejection is thus solely a matter of scale 
and feature shape, not intensity. Thus by deciding whether or not image points 
belong to a particular feature on the basis of feature asymmetry means that the 
technique is independent ofthe intensity of the image points. Of course the intensity 
can be noted and used in later processes, such as tracking as explained below. 

Figure 8 illustrates the results of processing an echocardiogram ofthe left ventricle 
with the log-Gabor functions of Figure 3 at a variety of scales from 16 pixels through 
to 72 pixels. The white overlay corresponds to detected step edges and marks those 
pixels with feature asymmetry values greater than the threshold 0.9. It can be seen 
that at the finest scale, for instance 16 pixels, most ofthe detected values correspond 
to the speckle texture in the ultrasound image. However, as the scale is increased it is 
the features of interest in the image which are detected, such as ventricular walls, for 
instance at 48 pixels scale. Thus for this image using a scale between 32 and 56 
pixels would give satisfectory results. The actual scale to be used in any situation 
depends on the image (for instance the quality of the image acquisition and the 
im a ging depth setting used by the sonographer), but the user can be given control of 
the scale giving the ability to tune or focus the detection to obtain good delineation of 
the desired features of interest. 

Figure 6 is an example that illustrates schematically the image analysis processing of 
this embodiment. The data 1 (which can be an image or sequence of images) is input 



wo 02/43004 



PCT/GBOl/05094 



'13- 

and Fourier transfonned at 3 before being convolved with the series of oriented 
quadrature filters 5. The output of the qxiadrature filters 5 is then used to compute the 
feature shape measure, which in this embodiment is the feature asymmetry (which is a 
maximum for step edges). The feature asymmetry is derived firom the even and odd 
5 filter outputs 7 using the formula above. The feature asjonmetry is compared to a 

threshold (for instance 0.9) which is the only parameter used to accept image pixels as 
feature candidates via AND gates 1 1 for subsequent processing. For processing in 
two spatial dimensions three filters are required, for processing in 3D six filters are 
needed as shown in Figures 7(A) and (B). In Figure 7(A) the orientations are directed 

10 to the vertices of a regular , hexagon. In Figure 7(B) the orientations are in the 

directions of the vertices of a regular icosahedron. For each of these candidate pixels 
an intensity-based measure, in this case the local amplitude as defined above, is 
computed at 9. Thus for each pixel accepted according to the feature asymmetry 
threshold, the local amplitude and feature asymmetry value from each of the oriented 

15 filters is recorded. 

It will be appreciated that the value of the local amplitude firom any of the filters 
depends on the relative orientation of the filter with respect to the feature. The 
maximum response is obtained when the filter is incident normal to the direction of 

20 the feature. Conversely, the minimum response is obtained when the filter is oriented 
parallel to the feature. In between, the local amplitude response decays smoothly and 
monotonically. This can be understood intuitively because the intensity profile along 
a line paralld to a step edge is constant, whereas the intensity profile along a line 
normal to the step edge gives the fiiU step edge change. Thus the interpolator 13 

25 det^mines the orientation of the feature from the local amplitude responses of the 
oriented filters 5. There are a variety of ways in which this can be done, such as a 
linear combination of the signals . from the different orientations or by considering the 
outputs of the filters as vectors havii^ a direction corresponding to the filter 
orientation and a length corresponding to the local amplitude, and fitting an eUipse to 

30 the end points of the vectors, the major axis of the fitted ellipse corresponding to the 
normal to the feature direction. The direction of the normal does not, of course, 
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distingmsh between positive and negative edges, but these are distingmshable by the 
sign of the odd filter output. 

A non-maximum suppression algorithm 14 can be run on the resiilts to narrow the 
5 detected edges to approximately pixel wide contour segments. Further, an intensity- 
related noise rejection, and preferably thresholding and cludn following, stage 16 is 
used to avoid streaking and to enhance the continuity of the contour features. The 
result is the output at 1 8 of a series of connected points which are detected as 
belonging to an edge feature, with a probability proportional to their local amplitude 
10 and feature asymmetry scores, and labelled with an estimate of the feature orientation. 
In the case of sequences of images, an additional estimate of the feature plane with 
respect to the temporal axis gives an estimate of feature image velocity through the 
sequence. 

15 ■ Thus having detected pixels which belong to a particular type of image feature (i.e. a 
feature defined by a particular shape in the intensity profile), and having labelled 
those pixels with the local amplitude, feature asymmetry, orientation and optionally 
' velocity, the system has provided a particularly rich description of the properties of 
the detected features. This rich description can be used to advantage in fi^rther 

20 processing of the irtiage, particularly when compared with techniques which just mark 
a pixel as belonging to a detected feature or not An example of such further 
processing will be described below, with reference to Figures 9 and 10. 

In the particular application considered here, to echocardiographic images of the left 
25 ventricle, the feature being detected is the left ventricular waU (endocardium) which 
is a closed contour having a shape as illustrated at 20 in Figure 9. As the heart beats 
the ventricle contracts as the ventricular wall moves and the contour detected in a 
sequence of images of the beating heart should follow the movement of the 
endocardium. While ideally it would be possible simply to detect the image features 
30 corresponding to the endocardium in each fi-ame of the sequence, in practice noise 

and artifacts tend to produce a disconnected representation of the cardiac boundaries. 
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and a number of spurious responses unrelated to the cardiac waUs. Tracking the 
feature through a sequence of images can interpolate the features detected in the 
individual frames to produce a continuous contour which follows, or tracks, the 
feature in the image. 

Various tracking techniques have been proposed, but tracking can be improved by 
using the rich description of the features provided by the detection technique above. 
This is illustrated in Figure 9. Firstly, taking a contour 20 inherited from an adjacent 
time frame, a search slab 22 is projected inwardly and outwardly from each point on 
the contour. Figure 9 illustrates an example in which these slabs are three pixels wide 
and twelve pixels long. Then each tangential row of three pixels within the slab is 
processed separately to find the strongest, local amplitude response of the three 
pixels, that agrees most closely with the orientation of the contour point and has the 
highest feature asymmetry score. This is repeated for each row of three pixels in the 
slab. In this embodiment one of the three pixels was selected as a possible candidate 
if it has the strongest local amplitude response amongst the three and has a feature 
asymmetry score greater than 0.8 and if the difference in orientation between the 
contour normal and the feature normal is less than 5 degrees.. Thus the candidates are 
selected making use of the labels attached during the feature detection process. It 
should be noted that the orientation measure takes into consideration the direction of 
the step edge (whether it is a positive or negative step). 

It should be noted that in the case of processing an image in two spatial dimmsions 
plus time, the velocity estimate for the points can be used to restrict the search space 
along the normal. Thus the length of the search slab is selected according to the 
velocity of the feature. 

The result of this is a onerpixel-wide slab 24 normal to the contour which consists of 
that pixel (if any) from each of the tangential rows in the search slab which satisfies 
the search criteria. 
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The next step is to draw the new contour using the results of the one-pixel wide 
normals 24. In this embodiment this is achieved by a dynamic programming 
technique which finds the optimal path linking a single point &om each of the 
normals 24, the optimal path being defined with reference to a cost function. The cost 
5 . function defines each transition between candidate points fit>m one normal 24 to the 
next as having a cost associated with it. Thus, referring to. Figure 10, the normals 24 
are stacked next to each other and for a given pixel k in the ith row, the cost 
associated with a transition to the corresponding pixel in the next row or each of its 
neighbours is calculated. 

10 

The cost function is designed to blend the local amplitude values and orientation 
information and it includes partial costs such as:- a large penalty cost if an entry on 
the normal 24 has been left blank because no feature has been found in the image; the 
cost related to the local amplitude measurement, the distance away from the previous 
15 contour point and a smoothing term. The optimal path is the one that links one point 
from each normal such that the sum of all costs is minimal. A least sqxiares 
transformation is then used to fit the contour inherited from the adjacent frame to the 
points selected by the dynamic programming routine. 

20 Of course other techniques for tracking the contour can be used, but the aim is to 
improve the tracking process using the labels attached during the feature detection 
process. 
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CLAIMS 

1 . A method of detecting features-of-interest in an image, comprising the steps 
of: 

making an intensity independent comparison of the image intensity profile 
^th a shape model for predefined features of interest; and 

storing for the profile a label comprisiag the level of similarity to the shape 

model. 

2. A method according to claim 1 wherein the label further comprises a measure 
of the orientation of the feature. 

3. A method according to claim 1 or 2 wherein the predefined features of 
interest comprise positive and negative step edges in the intensity profile, 

4. A method according to claim 1 , 2 or 3, further comprising performing the 
comparison on an intensity profile taken across successive frames of an image 
sequence to derive a measure of the velocity of image features moving through the 
sequence, and wherein said label ftirther comprises said measure of the velocity. 

5. A method according to claim 1, 2, 3, or 4, wherein the comparison comprises 
the steps of: 

convolving the intensity profile with a qiiadrature filter to derive the phase and 
amplitude of a pair of odd and even components of the profile in the spatial or spatio 
temporal frequency domains; 

deriving from the difference between the components a measure of one of the 
feature symmetry or asymmetry along the profile; 

deriving from the components a measure of the local amplitude along the 

profile; 

defining areas of the image satisfying a predefined constraint on the value of 
said one of the feature symmetry or asymmetry as satisfying the shape model and thus 
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belonging to an image feature; and 

vdierein the label comprises said measures. 

6. A method according to claim 5 wherein said quadrature filters are log-Gabor 
5 filters. 

7. A method according to claim 5 or 6 wherein said predefined constraint on said 
one of the feature synametry or asymmetry is that its magnitude is greater than or 
equal to a predefined threshold. 

10 

8. A method according to claim 5, 6, or 7, comprising the step of convolving the 
intensity profile of the image with quadrature filters oriented in different directions 
across the image to detect image features oriented differently in the image, and 
deriving from the relative responses of the filters a measure of the feature orientation. 

15 

9. • A method according to any one of the preceding claims, further comprising 
tracking detected image features through successive frames of an image sequence by 
comparing the labels of areas in a succeeding frame with the labels of areas in the 
preceding firame over a predetermined search region to define as relating to the same 

20 feature those areas in each fi:ame whose labels satisfy predefined search criteria. 

10. A method according to claim 9, wherein the predetemiined search region is 
defined as a region within a predetenniaed distance of the position, of the detected 
image feature in an adjacent frame in the sequence. 

25 

11. A method according to claim 10, wherein the predetemiined distance is set in 
dependence upon the velocity of the detected image feature. 

12. A method according to claim 9, 10 or 1 1 wherein the predefined search 

30 criteria comprise a' condition on the valiie of the feature symmetry or asymmetry of 
the areas in the search region. 



wo 02/43004 



PCT/GBOl/05094 



-19- 

13. A method according to claim 12 wherein the condition on the value of the 
feature symmetry or symmetry or asymmetry is that it is greater than a predeteraxined 
amoimt. 

5 14. A method according to claim 12 or 13 wherein the predefined search criteria 
comprise the condition that the difference in the orientation of the detected features in 
the two frames is less than a predetermined amount 

15. A method according to any one of the preceding claims wherein the intensity 
, 10 independent comparison with the shape model is performed at a plurality of scales. 

16. A method according to any one of the preceding claims wherein the 
comparison is performed on an intensity profile taken through a plurality of images of , 
different depths of a volume being imaged. 

15 

17. A method according to any one of the preceding claims wherein the image is 
an ultrasound or magnetic resonance or x-ray or standard camera image. 

18. A method according to any one of the preceding claims wherein the image is 
20 an ultrasound contrast image, an ultrasound Doppler image or an ultrasound harmonic 

image. 

19. A method au:cording to any one of. the preceding claims wherein the image is 
a medical image. 

25 

20. A method according to any one of the preceding claims wherein the image is a 
remote sensing or non-invasive image. 

21 . A method according to any one of claims 1 to 16 wherein the image is an 
30 echocardiogram. 
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22. A method according to claims 20 wherein the feature is a ventricular or atrial 
wall. 

23. A computer program comprising program code means for executing the 
5 method of any one of the preceding claims. 



24. An image processing system comprising image processing means adapted to 
perform the method of any one of the preceding claims. 
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Fig.10. 
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